# Replication materials for 
# Hoffmann, R., Muttarak, R., Peisker, J., Stanig, P.: 
# Climate change experiences raise environmental concerns and promote green voting

# set working directory
setwd("")
# this directory should contain  
# - eb_climate.RData -> environmental concern data set based on the Eurobarometer
# - vote_climate.RData -> Green voting data set of European Parliament elections
# - functions.R -> helper functions
# - conley.R -> correction of standard errors
# - iterate-obs-function.R -> helper functions for conley.R
# - cpp-functions.cpp -> helper function for conley.R
# - Tables folder 
# - Figures folder

# please run the code in the order in which it is provided

#### load packages and additional functions ####
library(tidyverse)
library(data.table)
library(lfe)
library(geosphere)
library(Rcpp)
library(RcppArmadillo)
library(stargazer)
library(cowplot)
library(countrycode)
library(xtable)
library(scico)
library(lmtest)
library(biostat3)
library(car)
library(parallel)

# citations
citation("tidyverse")
citation("data.table")
citation("lfe")
citation("geosphere")
citation("Rcpp")
citation("RcppArmadillo")
citation("stargazer")
citation("cowplot")
citation("countrycode")
citation("xtable")
citation("scico")
citation("lmtest")
citation("biostat3")
citation("car")
citation("parallel")

# load helper functions
source("functions.R")

# load functions for correction of standard errors by Darin Christensen
# Available at: https://github.com/darinchristensen/conley-se
# based on 
# Hsiang (2010); Temperatures and cyclones strongly associated with economic production in the Caribbean and Central America, 
# Proceedings of the National Academy of Sciences , Vol. 107, No. 35 
source("conley.R")

# set ggplot theme
theme_set(theme_bw())

print(sessionInfo(), locale = FALSE)
# R version 4.1.2 (2021-11-01)
# Platform: x86_64-pc-linux-gnu (64-bit)
# Running under: Arch Linux
# 
# Matrix products: default
# BLAS:   /usr/lib/libblas.so.3.10.0
# LAPACK: /usr/lib/liblapack.so.3.10.0
# 
# attached base packages:
# [1] parallel  stats     graphics  grDevices utils     datasets  methods   base     
# 
# other attached packages:
# [1] biostat3_0.1.5           MASS_7.3-54              survival_3.2-13          magrittr_2.0.1           car_3.0-10              
# [6] carData_3.0-4            lmtest_0.9-38            zoo_1.8-9                scico_1.2.0              xtable_1.8-4            
# [11] countrycode_1.2.0        cowplot_1.1.1            stargazer_5.2.2          RcppArmadillo_0.10.5.0.0 Rcpp_1.0.6              
# [16] geosphere_1.5-10         lfe_2.8-6                Matrix_1.3-4             data.table_1.14.0        forcats_0.5.1           
# [21] stringr_1.4.0            dplyr_1.0.6              purrr_0.3.4              readr_1.4.0              tidyr_1.1.3             
# [26] tibble_3.1.2             ggplot2_3.3.3            tidyverse_1.3.1         
# 
# loaded via a namespace (and not attached):
# [1] httr_1.4.2       jsonlite_1.7.2   splines_4.1.2    modelr_0.1.8     Formula_1.2-4    assertthat_0.2.1 sp_1.4-5        
# [8] cellranger_1.1.0 pillar_1.6.1     backports_1.2.1  lattice_0.20-45  glue_1.4.2       digest_0.6.27    rvest_1.0.0     
# [15] muhaz_1.2.6.4    colorspace_2.0-1 sandwich_3.0-1   survey_4.1-1     pkgconfig_2.0.3  broom_0.7.6      haven_2.4.1     
# [22] scales_1.1.1     openxlsx_4.2.3   rio_0.5.26       farver_2.1.0     generics_0.1.0   ellipsis_0.3.2   withr_2.4.2     
# [29] cli_2.5.0        crayon_1.4.1     readxl_1.3.1     fs_1.5.0         fansi_0.5.0      xml2_1.3.2       foreign_0.8-81  
# [36] tools_4.1.2      hms_1.1.0        mitools_2.4      lifecycle_1.0.0  munsell_0.5.0    reprex_2.0.0     zip_2.1.1       
# [43] compiler_4.1.2   rlang_0.4.11     grid_4.1.2       rstudioapi_0.13  labeling_0.4.2   gtable_0.3.0     abind_1.4-5     
# [50] DBI_1.1.1        curl_4.3.1       R6_2.5.0         lubridate_1.7.10 utf8_1.2.1       stringi_1.6.2    vctrs_0.3.8     
# [57] dbplyr_2.1.1     tidyselect_1.1.1

#### load data ####
# load voting data
load("vote_climate.RData")
# load environmental concern data
load("eb_climate.RData")

#### set parameters ####
# dependent variables
env_y <- "env_concern"
gv_y <- "green_vote"
# independent variables
xv <- c(
  "Temperature anomaly" = "tmean_nf_12m",
  "Heat episode (temp.)" = "tmean_xf_pos05_12m",
  "Heat episode (UTCI)" = "utci_mean_xf_pos05_12m",
  "Dry spell" = "spei3_neg_12m"
)
# fixed effects variables
env_fv <- c("unit", "year_cat", "season")
gv_fv <- c("unit", "year")
# cluster variables
env_cv <- c("cluster")
gv_cv <- c("cluster")
# auto-regressive terms
env_ar <- c("Concern (lag 1)" = "env_concern_l1")
gv_ar <- c("Green vote (lag 1)" = "green_vote_l1")
# time-varying control variables 
controls <- c(
  "GDP per capita" = "gdp_ln_l6", 
  "Unemployment rate" = "unemp_lfs_l6", 
  "Agricultural share" = "gva_agriculture_share_l6"
  )
# disaster control variables
disasters <- c(
  "Flooding" = "flood_frac_12m",
  "Wildfire" = "burn_frac_12m"
)

# cutoff values for Conley correction of standard errors
# spatial distance in kilometers
env_spc <- 500
gv_spc <- 500
# temporal distance
env_tc <- 18 # in months
gv_tc <- 5 # in years
# set number of threads used, defaults to 4
# conley_threads <- 4

#### Table 1: Baseline effects of climate extremes on environmental concerns and Green voting ####
tab1 <- c(
  lapply(xv, function(i){ felm(
    f_felm(yvar = env_y, xvars = i, fevars = env_fv, clustervars = env_cv), 
    eb_climate, keepCX = TRUE, keepX = TRUE
  )}), 
  lapply(xv, function(i){ felm(
    f_felm(yvar = gv_y, xvars = i, fevars = gv_fv, clustervars = gv_cv), 
    vote_climate, keepCX = TRUE, keepX = TRUE
  )})
)
tab1_cs <- c(
  lapply(tab1[1:length(xv)], conley_concern),
  lapply(tab1[(1:length(xv))+length(xv)], conley_vote)
)
export_html(lapply(tab1_cs, standardise), out = "Tables/table1.html")

#### Table S1: Data summary for the concern dataset ####
paste0("Number of countries in concern dataset: ", length(unique(eb_climate$nuts0)))
paste0("Number of regions in concern dataset: ", length(unique(eb_climate$unit)))

table_s1 <- 
  eb_climate[
    !is.na(env_concern), 
    .(
      Country = unique(countrycode(nuts0, "eurostat", "country.name")),
      N = .N,
      Region = unique(region),
      `# NUTS` = length(unique(unit)),
      `# cold NUTS` = length(unique(unit[climate == "Cold"])),
      `# temperate NUTS` = length(unique(unit[climate == "Temperate"])),
      `# hot NUTS` = length(unique(unit[climate == "Hot"])),
      `Series start` = min(year),
      `# Eurobarometer` = as.integer(max(table(unit)))
    ),
    by = .(nuts0)
  ] %>% 
  .[order(Region, Country), -1] 

print(
  xtable::xtable(table_s1, digits = 1),
  file = "Tables/table_s1.html",
  include.rownames = FALSE,
  type = "html"
)

#### Table S2: Data summary for the voting dataset ####
paste0("Number of countries in election dataset: ", length(unique(vote_climate$nuts0)))
paste0("Number of regions in election dataset: ", length(unique(vote_climate$unit)))

table_s2 <- 
  vote_climate[
    !is.na(green_vote), 
    .(
      Country = unique(countrycode(nuts0, "eurostat", "country.name")),
      N = .N,
      Region = unique(region),
      `# NUTS` = length(unique(unit)),
      `# cold NUTS` = length(unique(unit[climate == "Cold"])),
      `# temperate NUTS` = length(unique(unit[climate == "Temperate"])),
      `# hot NUTS` = length(unique(unit[climate == "Hot"])),
      `Series start` = min(year),
      `# elections` = as.integer(max(table(unit)))
    ),
    by = .(nuts0)
  ] %>% 
  .[order(Region, Country), -1] 

print(
  xtable::xtable(table_s2, digits = 1),
  file = "Tables/table_s2.html",
  include.rownames = FALSE,
  type = "html"
)

#### Table S3: Tests for serial correlation of residuals ####
# get residuals from table 1
tab1_res_concern <- c(
  lapply(seq_along(xv), function(i){ return(
    eb_climate %>% 
      select(unit, t, env_dt1, env_dt2, env_dt3, env_dt4) %>% 
      bind_cols(res = tab1[[i]]$residuals) %>% 
      group_by(unit) %>% 
      mutate(res_l1 = shift(res, 1), res_l2 = shift(res, 2), res_l3 = shift(res, 3), res_l4 = shift(res, 4), res_l5 = shift(res, 5)) %>% 
      ungroup()
  )}))
tab1_res_vote <- c(
  lapply(seq_along(xv), function(i){ return(
    vote_climate %>% 
      select(unit, year) %>% 
      bind_cols(res = tab1[[i+length(xv)]]$residuals) %>% 
      group_by(unit) %>% 
      mutate(res_l1 = shift(res, 1), res_l2 = shift(res, 2), res_l3 = shift(res, 3)) %>% 
      ungroup()
  )}))

# regress residuals on lagged values
res_reg <- c(
  lapply(tab1_res_concern, function(x){
    lm(res ~ 
         res_l1 + res_l1:env_dt1 + 
         res_l2 + res_l2:env_dt2 + 
         res_l3 + res_l3:env_dt3 + 
         res_l4 + res_l4:env_dt4, 
       data = x)
  }),
  lapply(tab1_res_vote, function(x){
    lm(res ~ res_l1 + res_l2, data = x)
  })
)
export_html(
  res_reg, out = "Tables/table_s3.html",
  cov_labels = c(
    "Residual (lag 1)", "Residual (lag 2)", "Residual (lag 3)","Residual (lag 4)",
    "Residual (lag 1) × dt1", "Residual (lag 2) × dt2", "Residual (lag 3) × dt3","Residual (lag 4) × dt4"
  ),
  dep_labels = "Residual",
  col_labels = c("Environmental concern", "Green vote share"),
  add_lines = NULL,
  omit_stats = c("ser", "adj.rsq", "f")
)

#### Table S4: Permutation tests for Moran’s I statistic of the residuals ####
# available on request

#### Table S5: Summary statistics of the concern dataset ####
sum_var_labels <- c(
  "env_concern" = "Environmental concern",
  "tmean_nf_12m" = "Temperature anomaly",
  "tmean_nf_pos_12m" = "Temperature anomaly (positive)",
  "tmean_nf_neg_12m" = "Temperature anomaly (negative)",
  "tmean_xf_pos05_12m" = "Heat episode (temp.)",
  "tmean_xf_neg05_12m" = "Cold episode (temp.)",
  "utci_mean_xf_pos05_12m" = "Heat episode (UTCI)",
  "utci_mean_xf_neg05_12m" = "Cold episode (UTCI)",
  "spei3_neg_12m" = "Dry spell (SPEI3)",
  "spei3_pos_12m" = "Wet spell (SPEI3)",
  "gdp_ln_within" = "Log GDP per capita (within)",
  "gdp_ln_mean" = "Log GDP per capita (between)",
  "agri_frac" = "Agricultural share of region (between)",
  "irrigated_mean" = "Share of agr. area with irrigation (between)",
  "share_edu3_mean" = "Pop. share with at least secondary education (between)",
  "pop_urban" = "Urban region (between)",
  "LegislativeFractionalization_mean" = "Legislative Fractionalization (between)",
  "CenterOfGravity_1" = "Incumbent’s ideology"
)
sum_vars <- names(sum_var_labels)
var_names <- 
  sum_var_labels %>% 
  data.frame %>% 
  rownames_to_column()
names(var_names) <- c("variable", "Variable")

# get variance of variables within units
reg_vars <- sum_vars[!grepl("between", sum_var_labels, fixed = TRUE)]
within_sd_concern <- 
  sapply(reg_vars, function(i){ 
    felm_tmp <- felm(  
      f_felm(yvar = "env_concern", xvars = i, fevars = env_fv), 
      eb_climate, keepCX = TRUE, keepX = TRUE)
    sd_tmp <- c(sd(felm_tmp$cX))
    return(sd_tmp)
  }) %>% 
  data.frame() %>% 
  rownames_to_column 
names(within_sd_concern) <- c("variable", "Within SD")

# calculate summary statistics
round_cols <- c("Min", "Mean", "Median", "Max", "SD", "Within SD")
table_s5 <- 
  eb_climate[, ..sum_vars] %>% 
  melt(measure.vars = 1:ncol(.)) %>% 
  .[
    , .(
      N = sum(!is.na(value)),
      Min = min(value, na.rm = TRUE),
      Mean = mean(value, na.rm = TRUE),
      Median = median(value, na.rm = TRUE),
      Max = max(value, na.rm = TRUE),
      SD = sd(value, na.rm = TRUE)
    ),
    by = .(variable)
  ] %>% 
  merge(within_sd_concern, by = "variable", all.x = TRUE) %>% 
  merge(var_names, by = "variable", all.x = TRUE) %>% 
  .[
    , c(round_cols) := lapply(.SD, function(x){round(x, 3)}),
    .SDcols = c(round_cols)
  ] %>% 
  .[match(sum_vars, variable), c(9, 2:8)]

print(
  xtable::xtable(table_s5, digits = 3), 
  file = "Tables/table_s5.html",
  include.rownames = FALSE,
  type = "html"
)

#### Table S6: Summary statistics of the voting dataset ####
sum_var_labels <- c(
  "green_vote" = "Green vote share",
  "tmean_nf_12m" = "Temperature anomaly",
  "tmean_nf_pos_12m" = "Temperature anomaly (positive)",
  "tmean_nf_neg_12m" = "Temperature anomaly (negative)",
  "tmean_xf_pos05_12m" = "Heat episode (temp.)",
  "tmean_xf_neg05_12m" = "Cold episode (temp.)",
  "utci_mean_xf_pos05_12m" = "Heat episode (UTCI)",
  "utci_mean_xf_neg05_12m" = "Cold episode (UTCI)",
  "spei3_neg_12m" = "Dry spell (SPEI3)",
  "spei3_pos_12m" = "Wet spell (SPEI3)",
  "gdp_ln_within" = "Log GDP per capita (within)",
  "gdp_ln_mean" = "Log GDP per capita (between)",
  "agri_frac" = "Agricultural share of region (between)",
  "irrigated_mean" = "Share of agr. area with irrigation (between)",
  "share_edu3_mean" = "Pop. share with at least secondary education (between)",
  "pop_urban" = "Urban region (between)",
  "LegislativeFractionalization_mean" = "Legislative Fractionalization (between)",
  "CenterOfGravity_1" = "Incumbent’s ideology"
)
sum_vars <- names(sum_var_labels)
var_names <- 
  sum_var_labels %>% 
  data.frame %>% 
  rownames_to_column()
names(var_names) <- c("variable", "Variable")

# get variance of variables within units
reg_vars <- sum_vars[!grepl("between", sum_var_labels, fixed = TRUE)]
within_sd_vote <- 
  sapply(reg_vars, function(i){ 
    felm_tmp <- felm(  
      f_felm(yvar = "green_vote", xvars = i, fevars = gv_fv), 
      vote_climate, keepCX = TRUE, keepX = TRUE)
    sd_tmp <- c(sd(felm_tmp$cX))
    return(sd_tmp)
  }) %>% 
  data.frame() %>% 
  rownames_to_column 
names(within_sd_vote) <- c("variable", "Within SD")
# save(within_sd_vote, file = "Results/R files/within_sd_vote.RData")

# calculate summary statistics
round_cols <- c("Min", "Mean", "Median", "Max", "SD", "Within SD")
table_s6 <- 
  vote_climate[, ..sum_vars] %>% 
  melt(measure.vars = 1:ncol(.)) %>% 
  .[
    , .(
      N = sum(!is.na(value)),
      Min = min(value, na.rm = TRUE),
      Mean = mean(value, na.rm = TRUE),
      Median = median(value, na.rm = TRUE),
      Max = max(value, na.rm = TRUE),
      SD = sd(value, na.rm = TRUE)
    ),
    by = .(variable)
  ] %>% 
  merge(within_sd_vote, by = "variable", all.x = TRUE) %>% 
  merge(var_names, by = "variable", all.x = TRUE) %>% 
  .[
    , c(round_cols) := lapply(.SD, function(x){round(x, 3)}),
    .SDcols = c(round_cols)
  ] %>% 
  .[match(sum_vars, variable), c(9, 2:8)]

print(
  xtable::xtable(table_s6, digits = 3), 
  file = "Tables/table_s6.html",
  include.rownames = FALSE,
  type = "html"
)

#### Table S7: Baseline models with clustered standard errors ####
# see Table 1
export_html(
  lapply(tab1, standardise), 
  out = "Tables/table_s7.html",
  add_lines = list(
    c("Unit fixed effects", rep("×", length(xv)*2)),
    c("Period fixed effects", rep("×", length(xv)*2)),
    c("Season fixed effects", rep("×", length(xv)*2)),
    c("Cluster", rep(toupper(env_cv), length(xv)), rep(toupper(gv_cv), length(xv))),
    c("Observations", sapply(tab1_cs, get_n)),
    c("Overall R²", sapply(tab1_cs, get_r2_overall)),
    c("Within R²", sapply(tab1_cs, get_r2_within))
  ))

#### Table S8: Baseline models controlling for lagged dependent variables ####
table_s8 <- c(
  lapply(xv, function(i){ felm(
    f_felm(yvar = env_y, xvars = c(i, env_ar), fevars = env_fv, clustervars = env_cv), 
    eb_climate, keepCX = TRUE, keepX = TRUE
  )}), 
  lapply(xv, function(i){ felm(
    f_felm(yvar = gv_y, xvars = c(i, gv_ar), fevars = gv_fv, clustervars = gv_cv), 
    vote_climate, keepCX = TRUE, keepX = TRUE
  )}) 
)
table_s8_cs <- c(
  lapply(table_s8[1:length(xv)], conley_concern),
  lapply(table_s8[(1:length(xv))+length(xv)], conley_vote)
)
export_html(lapply(table_s8_cs, standardise), out = "Tables/table_s8.html", cov_labels = c(names(xv), names(env_ar), names(gv_ar)))

#### Table S9: Baseline models restricted to the time period of 2007–2019 ####
table_s9 <- c(
  lapply(xv, function(i){ felm(
    f_felm(yvar = env_y, xvars = i, fevars = env_fv, clustervars = env_cv), 
    eb_climate[year >= 2007, ], keepCX = TRUE, keepX = TRUE
  )}), 
  lapply(xv, function(i){ felm(
    f_felm(yvar = gv_y, xvars = i, fevars = gv_fv, clustervars = gv_cv), 
    vote_climate[year >= 2007, ], keepCX = TRUE, keepX = TRUE
  )}) 
)
table_s9_cs <- c(
  lapply(table_s9[1:length(xv)], conley_concern, data = eb_climate[year >= 2007, ]), 
  lapply(table_s9[(1:length(xv))+length(xv)], conley_vote, data = vote_climate[year >= 2007, ])
)
export_html(lapply(table_s9_cs, standardise), out = "Tables/table_s9.html")

#### Table S10: Robustness test: Baseline models restricted to the time period of 2011–2019 ####
table_s10 <- c(
  lapply(xv, function(i){ felm(
    f_felm(yvar = env_y, xvars = i, fevars = env_fv, clustervars = env_cv), 
    eb_climate[year >= 2011, ], keepCX = TRUE, keepX = TRUE
  )}), 
  lapply(xv, function(i){ felm(
    f_felm(yvar = gv_y, xvars = i, fevars = gv_fv, clustervars = gv_cv), 
    vote_climate[year >= 2011, ], keepCX = TRUE, keepX = TRUE
  )}) 
)
table_s10_cs <- c(
  lapply(table_s10[1:length(xv)], conley_concern, data = eb_climate[year >= 2011, ]), 
  lapply(table_s10[(1:length(xv))+length(xv)], conley_vote, data = vote_climate[year >= 2011, ])
)
export_html(lapply(table_s10_cs, standardise), out = "Tables/table_s10.html")

#### Table S11: Baseline models with climate change concern as outcome ####
# available on request (requires additional data)

#### Table S12: Baseline models with time-varying controls ####
table_s12 <- c(
  lapply(xv, function(i){ felm(
    f_felm(yvar = env_y, xvars = c(i, controls), fevars = env_fv, clustervars = env_cv), 
    eb_climate, keepCX = TRUE, keepX = TRUE
  )}), 
  lapply(xv, function(i){ felm(
    f_felm(yvar = gv_y, xvars = c(i, controls), fevars = gv_fv, clustervars = gv_cv), 
    vote_climate, keepCX = TRUE, keepX = TRUE
  )}) 
)
table_s12_cs <- c(
  lapply(table_s12[1:length(xv)], conley_concern),
  lapply(table_s12[(1:length(xv))+length(xv)], conley_vote)
)
export_html(lapply(table_s12_cs, standardise), out = "Tables/table_s12.html", cov_labels = c(names(xv), names(controls)))

#### Table S13: Baseline models using 1961–1990 as climate reference period ####
xv <- c(
  "tmean_nfb_12m",
  "tmean_nfb_pos_12m + tmean_nfb_neg_12m",
  "tmean_xf_pos05b_12m + tmean_xf_neg05b_12m",
  "spei61_neg_12m + spei61_pos_12m"
)
names_xv <- c(
  "Temperature anomaly",
  "Temperature anomaly (+)", "Temperature anomaly (–)",
  "Heat episode (temp.)", "Cold episode (temp.)",
  "Dry spell", "Wet spell"
)
table_s13 <- c(
  lapply(xv, function(i){ felm(
    f_felm(yvar = env_y, xvars = i, fevars = env_fv, clustervars = env_cv), 
    eb_climate, keepCX = TRUE, keepX = TRUE
  )}), 
  lapply(xv, function(i){ felm(
    f_felm(yvar = gv_y, xvars = i, fevars = gv_fv, clustervars = gv_cv), 
    vote_climate, keepCX = TRUE, keepX = TRUE
  )}) 
)
table_s13_cs <- c(
  lapply(table_s13[1:length(xv)], conley_concern),
  lapply(table_s13[(1:length(xv))+length(xv)], conley_vote)
)
export_html(lapply(table_s13_cs, standardise), out = "Tables/table_s13.html", cov_labels = names_xv)

#### Table S14: Impacts of heat and cold episodes using different thresholds to define climate extremes ####
xv <- c(
  "tmean_xf_pos025_12m + tmean_xf_neg025_12m",
  "tmean_xf_pos05_12m + tmean_xf_neg05_12m",
  "tmean_xf_pos10_12m + tmean_xf_neg10_12m",
  "utci_mean_xf_pos025_12m + utci_mean_xf_neg025_12m",
  "utci_mean_xf_pos05_12m + utci_mean_xf_neg05_12m",
  "utci_mean_xf_pos10_12m + utci_mean_xf_neg10_12m"
)
names_xv <- c(
  "Heat episode (2.5% temp.)", "Cold episode (2.5% temp.)",
  "Heat episode (5% temp.)", "Cold episode (5% temp.)",
  "Heat episode (10% temp.)", "Cold episode (10% temp.)",
  "Heat episode (2.5% UTCI)", "Cold episode (2.5% UTCI)",
  "Heat episode (5% UTCI)", "Cold episode (5% UTCI)",
  "Heat episode (10% UTCI)", "Cold episode (10% UTCI)"
)
table_s14 <- c(
  lapply(xv, function(i){ felm(
    f_felm(yvar = env_y, xvars = i, fevars = env_fv, clustervars = env_cv), 
    eb_climate, keepCX = TRUE, keepX = TRUE
  )}), 
  lapply(xv, function(i){ felm(
    f_felm(yvar = gv_y, xvars = i, fevars = gv_fv, clustervars = gv_cv), 
    vote_climate, keepCX = TRUE, keepX = TRUE
  )}) 
)
table_s14_cs <- c(
  lapply(table_s14[1:length(xv)], conley_concern),
  lapply(table_s14[(1:length(xv))+length(xv)], conley_vote)
)
export_html(lapply(table_s14_cs, standardise), out = "Tables/table_s14.html", cov_labels = names_xv)

#### Table S15: Baseline concern models with post-stratification weighting ####
xv <- c(
  "Temperature anomaly" = "tmean_nf_12m",
  "Heat episode (temp.)" = "tmean_xf_pos05_12m",
  "Heat episode (UTCI)" = "utci_mean_xf_pos05_12m",
  "Dry spell" = "spei3_neg_12m"
)
table_s15 <- c(
  lapply(xv, function(i){ felm(
    f_felm(yvar = "env_concern_w", xvars = i, fevars = env_fv, clustervars = env_cv), 
    eb_climate, keepCX = TRUE, keepX = TRUE
  )})
)
table_s15_cs <- c(
  lapply(tab1[1:length(xv)], conley_concern),
  lapply(table_s15[1:length(xv)], conley_concern)
)
export_html(lapply(table_s15_cs, standardise), out = "Tables/table_s15.html",
            dep_label_include = FALSE,
            col_labels = c("Env. concern (unweighted)", "Env. concern (weighted)"),
            col_sep = rep(length(xv), 2),
            add_lines = list(
              c("Unit fixed effects", rep("×", length(xv)*2)),
              c("Period fixed effects", rep("×", length(xv)*2)),
              c("Season fixed effects", rep("×", length(xv)*2)),
              c("Spatial cutoff (km)", rep(env_spc, length(xv)*2)),
              c("Temporal cutoff (years)", rep(env_tc/12, length(xv)*2)),
              c("Observations", sapply(table_s15_cs, get_n)),
              c("Overall R²", sapply(table_s15_cs, get_r2_overall)),
              c("Within R²", sapply(table_s15_cs, get_r2_within))
            ))

#### Table S16: Baseline models controlling for wildfire events ####
table_s16 <- c(
  lapply(xv, function(i){ felm(
    f_felm(yvar = env_y, xvars = c(i, disasters[2]), fevars = env_fv, clustervars = env_cv), 
    eb_climate, keepCX = TRUE, keepX = TRUE
  )}),
  lapply(xv, function(i){ felm(
    f_felm(yvar = gv_y, xvars = c(i, disasters[2]), fevars = gv_fv, clustervars = gv_cv), 
    vote_climate, keepCX = TRUE, keepX = TRUE
  )}) 
)
table_s16_cs <- c(
  lapply(table_s16[1:length(xv)], conley_concern),
  lapply(table_s16[(1:length(xv))+length(xv)], conley_vote)
)
export_html(lapply(table_s16_cs, standardise), out = "Tables/table_s16.html", cov_labels = c(names(xv), names(disasters[2])))

#### Table S17: Baseline models controlling for flooding events ####
table_s17 <- c(
  lapply(xv, function(i){ felm(
    f_felm(yvar = env_y, xvars = c(i, disasters[1]), fevars = env_fv, clustervars = env_cv), 
    eb_climate, keepCX = TRUE, keepX = TRUE
  )}),
  lapply(xv, function(i){ felm(
    f_felm(yvar = gv_y, xvars = c(i, disasters[1]), fevars = gv_fv, clustervars = gv_cv), 
    vote_climate, keepCX = TRUE, keepX = TRUE
  )}) 
)
table_s17_cs <- c(
  lapply(table_s17[1:length(xv)], conley_concern),
  lapply(table_s17[(1:length(xv))+length(xv)], conley_vote)
)
export_html(lapply(table_s17_cs, standardise), out = "Tables/table_s17.html", cov_labels = c(names(xv), names(disasters[1])))

#### Table S18: Baseline models with unstandardized outcome and standardized independent variables ####
# see table 1
export_html(lapply(tab1_cs, standardise, stan_y = FALSE), out = "Tables/table_s18.html")

#### Table S19: Baseline models without standardization ####
# see table 1
export_html(tab1_cs, out = "Tables/table_s19.html")

#### Table S20: Baseline models standardized using the full variance of the outcomes and the residual variance of the regressors ####
# see table 1
export_html(lapply(tab1_cs, standardise, res_y = FALSE), out = "Tables/table_s20.html")

#### Table S21: Effects of positive and negative temperature extremes on concerns and voting ####
xv <- c(
  "tmean_nf_pos_12m + tmean_nf_neg_12m",
  "tmean_xf_pos05_12m + tmean_xf_neg05_12m",
  "utci_mean_xf_pos05_12m + utci_mean_xf_neg05_12m",
  "spei3_neg_12m + spei3_pos_12m"
)
names_xv <- c(
  "Temperature anomaly (+)", "Temperature anomaly (–)",
  "Heat episode (temp.)", "Cold episode (temp.)",
  "Heat episode (UTCI)", "Cold episode (UTCI)",
  "Dry spell", "Wet spell"
)
table_s21 <- c(
  lapply(xv, function(i){ felm(
    f_felm(yvar = env_y, xvars = i, fevars = env_fv, clustervars = env_cv), 
    eb_climate, keepCX = TRUE, keepX = TRUE
  )}), 
  lapply(xv, function(i){ felm(
    f_felm(yvar = gv_y, xvars = i, fevars = gv_fv, clustervars = gv_cv), 
    vote_climate, keepCX = TRUE, keepX = TRUE
  )})
)
table_s21_cs <- c(
  lapply(table_s21[1:length(xv)], conley_concern),
  lapply(table_s21[(1:length(xv))+length(xv)], conley_vote)
)
export_html(lapply(table_s21_cs, standardise), out = "Tables/table_s21.html", cov_labels = names_xv)

#### Table S22: Effects of temperature anomalies on concerns and voting with different time lags ####
xv <- c(
  "tmean_nf_1m",
  "tmean_nf_6m",
  "tmean_nf_12m",
  "tmean_nf_24m",
  "tmean_nf_36m",
  "tmean_nf_48m"
)
names_xv <- c(
  "Temperature anomaly (1 month)",
  "Temperature anomaly (6 months)",
  "Temperature anomaly (12 months)",
  "Temperature anomaly (24 months)",
  "Temperature anomaly (36 months)",
  "Temperature anomaly (48 months)"
)
table_s22 <- c(
  lapply(xv, function(i){ felm(
    f_felm(yvar = env_y, xvars = i, fevars = env_fv, clustervars = env_cv), 
    eb_climate, keepCX = TRUE, keepX = TRUE
  )}), 
  lapply(xv, function(i){ felm(
    f_felm(yvar = gv_y, xvars = i, fevars = gv_fv, clustervars = gv_cv), 
    vote_climate, keepCX = TRUE, keepX = TRUE
  )}) 
)
table_s22_cs <- c(
  lapply(table_s22[1:length(xv)], conley_concern),
  lapply(table_s22[(1:length(xv))+length(xv)], conley_vote)
)
export_html(lapply(table_s22_cs, standardise), out = "Tables/table_s22.html", cov_labels = names_xv)

#### Table S23: Effects of heat episodes (temperature-based) on concerns and voting with different time lags ####
xv <- c(
  "tmean_xf_pos05_1m",
  "tmean_xf_pos05_6m",
  "tmean_xf_pos05_12m",
  "tmean_xf_pos05_24m",
  "tmean_xf_pos05_36m",
  "tmean_xf_pos05_48m"
)
names_xv <- c(
  "Heat episode (1 month, temp.)",
  "Heat episode (6 months, temp.)",
  "Heat episode (12 months, temp.)",
  "Heat episode (24 months, temp.)",
  "Heat episode (36 months, temp.)",
  "Heat episode (48 months, temp.)"
)
table_s23 <- c(
  lapply(xv, function(i){ felm(
    f_felm(yvar = env_y, xvars = i, fevars = env_fv, clustervars = env_cv), 
    eb_climate, keepCX = TRUE, keepX = TRUE
  )}), 
  lapply(xv, function(i){ felm(
    f_felm(yvar = gv_y, xvars = i, fevars = gv_fv, clustervars = gv_cv), 
    vote_climate, keepCX = TRUE, keepX = TRUE
  )}) 
)
table_s23_cs <- c(
  lapply(table_s23[1:length(xv)], conley_concern),
  lapply(table_s23[(1:length(xv))+length(xv)], conley_vote)
)
export_html(lapply(table_s23_cs, standardise), out = "Tables/table_s23.html", cov_labels = names_xv)

#### Table S24: Effects of heat episodes (UTCI-based) on concerns and voting with different time lags ####
xv <- c(
  "utci_mean_xf_pos05_1m",
  "utci_mean_xf_pos05_6m",
  "utci_mean_xf_pos05_12m",
  "utci_mean_xf_pos05_24m",
  "utci_mean_xf_pos05_36m",
  "utci_mean_xf_pos05_48m"
)
names_xv <- c(
  "Heat episode (1 month, UTCI)",
  "Heat episode (6 months, UTCI)",
  "Heat episode (12 months, UTCI)",
  "Heat episode (24 months, UTCI)",
  "Heat episode (36 months, UTCI)",
  "Heat episode (48 months, UTCI)"
)
table_s24 <- c(
  lapply(xv, function(i){ felm(
    f_felm(yvar = env_y, xvars = i, fevars = env_fv, clustervars = env_cv), 
    eb_climate, keepCX = TRUE, keepX = TRUE
  )}), 
  lapply(xv, function(i){ felm(
    f_felm(yvar = gv_y, xvars = i, fevars = gv_fv, clustervars = gv_cv), 
    vote_climate, keepCX = TRUE, keepX = TRUE
  )}) 
)
table_s24_cs <- c(
  lapply(table_s24[1:length(xv)], conley_concern),
  lapply(table_s24[(1:length(xv))+length(xv)], conley_vote)
)
export_html(lapply(table_s24_cs, standardise), out = "Tables/table_s24.html", cov_labels = names_xv)

#### Table S25: Effects of dry spells on concerns and voting with different time lags ####
xv <- c(
  "spei3_neg_1m",
  "spei3_neg_6m",
  "spei3_neg_12m",
  "spei3_neg_24m",
  "spei3_neg_36m",
  "spei3_neg_48m"
)
names_xv <- c(
  "Dry spell (1 month, SPEI3)",
  "Dry spell (6 months, SPEI3)",
  "Dry spell (12 months, SPEI3)",
  "Dry spell (24 months, SPEI3)",
  "Dry spell (36 months, SPEI3)",
  "Dry spell (48 months, SPEI3)"
)
table_s25 <- c(
  lapply(xv, function(i){ felm(
    f_felm(yvar = env_y, xvars = i, fevars = env_fv, clustervars = env_cv), 
    eb_climate, keepCX = TRUE, keepX = TRUE
  )}), 
  lapply(xv, function(i){ felm(
    f_felm(yvar = gv_y, xvars = i, fevars = gv_fv, clustervars = gv_cv), 
    vote_climate, keepCX = TRUE, keepX = TRUE
  )}) 
)
table_s25_cs <- c(
  lapply(table_s25[1:length(xv)], conley_concern),
  lapply(table_s25[(1:length(xv))+length(xv)], conley_vote)
)
export_html(lapply(table_s25_cs, standardise), out = "Tables/table_s25.html", cov_labels = names_xv)

#### Table S26: Climatic impacts on concerns and voting by different climatic zones ####
xv <- c(
  "tmean_nf_12m + tmean_nf_12m:climate",
  "tmean_xf_pos05_12m + tmean_xf_pos05_12m:climate",
  "utci_mean_xf_pos05_12m + utci_mean_xf_pos05_12m:climate",
  "spei3_neg_12m + spei3_neg_12m:climate"
)
names_xv <- c(
  "Temperature anomaly", "Temperature anomaly × Hot", "Temperature anomaly × Temperate",
  "Heat episode (temp.)", "Heat episode (temp.) × Hot", "Heat episode (temp.) × Temperate",
  "Heat episode (UTCI)", "Heat episode (UTCI) × Hot", "Heat episode (UTCI) × Temperate",
  "Dry spell", "Dry spell × Hot", "Dry spell × Temperate"
)
table_s26 <- c(
  lapply(xv, function(i){ felm(
    f_felm(yvar = env_y, xvars = i, fevars = env_fv, clustervars = env_cv), 
    eb_climate, keepCX = TRUE, keepX = TRUE
  )}), 
  lapply(xv, function(i){ felm(
    f_felm(yvar = gv_y, xvars = i, fevars = gv_fv, clustervars = gv_cv),
    vote_climate, keepCX = TRUE, keepX = TRUE
  )})
)
table_s26_cs <- c(
  lapply(table_s26[1:length(xv)], conley_concern),
  lapply(table_s26[(1:length(xv))+length(xv)], conley_vote)
)
export_html(table_s26_cs, out = "Tables/table_s26.html", cov_labels = names_xv)

#### Figure 2: Dot-whisker plot of temperature effects by climate ####
lincom_concern <- 
  lapply(table_s26_cs[1:length(xv)], function(x){
    coef_names <- rownames(x$coefficients)
    # get SD of outcome by climate type 
    sd_y <- data.table(x$cY, climate = eb_climate$climate)[
      , .(sd_y = sd(env_concern)), by = .(climate)
    ]
    # get SD of covariates by climate type
    sd_x <- data.table(x$cX, climate = eb_climate$climate)[
      , lapply(.SD, sd),
      by = .(climate),
      .SDcols = c(coef_names[!grepl(":", coef_names)])
    ] %>% 
      pivot_longer(2:ncol(.), names_to = "term", values_to = "sd_x")
    # calculate linear combination for each climate
    effects <- bind_rows(
      biostat3::lincom(x, paste0(coef_names[1])),
      biostat3::lincom(x, paste0(c(coef_names[1], coef_names[2]), collapse = "+")),
      biostat3::lincom(x, paste0(c(coef_names[1], coef_names[3]), collapse = "+"))
    ) %>% 
      bind_rows() %>% rownames_to_column() %>% 
      separate(rowname, into = c("term", "climate"), sep = "\\+") %>% 
      mutate(
        climate = case_when(
          grepl("Hot", climate) ~ "Hot",
          grepl("Temperate", climate) ~ "Temperate",
          TRUE ~ "Cold"
        )
      )
    # standardise by climate
    effects_stan <- 
      effects %>% 
      merge(sd_y, by = "climate") %>% 
      merge(sd_x, by = c("climate", "term")) %>% 
      transmute(
        term, 
        Climate = factor(climate, levels = c("Hot", "Temperate", "Cold")),
        Estimate = Estimate * (sd_x / sd_y),
        conf.low = `2.5 %`  * (sd_x / sd_y),
        conf.high = `97.5 %`  * (sd_x / sd_y),
        statistic = Chisq, 
        p.value = `Pr(>Chisq)`
      )
    return(effects_stan)
  }) %>% 
  # combine different models into one data frame
  bind_rows() %>% 
  mutate(
    term = case_when(
      term %in% c("tmean_nf_12m") ~ "Temperature anomaly",
      term %in% c("tmean_nf_pos_12m") ~ "Temperature anomaly (+)",
      term %in% c("tmean_nf_neg_12m") ~ "Temperature anomaly (–)",
      term %in% c("tmean_xf_pos05_12m") ~ "Heat episode (temp.)",
      term %in% c("tmean_xf_neg05_12m") ~ "Cold episode (temp.)",
      term %in% c("utci_mean_xf_pos05_12m") ~ "Heat episode (UTCI)",
      term %in% c("utci_mean_xf_neg05_12m") ~ "Cold episode (UTCI)",
      term %in% c("spei3_neg_12m") ~ "Dry spell",
      term %in% c("spei3_pos_12m") ~ "Wet spell"
    ),
    dep_var = "Environmental concern"
  )

lincom_vote <- 
  lapply(table_s26_cs[(1:length(xv))+length(xv)], function(x){
    coef_names <- rownames(x$coefficients)
    # get SD of outcome by climate 
    sd_y <- data.table(x$cY, climate = vote_climate$climate)[
      , .(sd_y = sd(green_vote)), by = .(climate)
    ]
    # get SD of covariates by climate 
    sd_x <- data.table(x$cX, climate = vote_climate$climate)[
      , lapply(.SD, sd),
      by = .(climate),
      .SDcols = c(coef_names[!grepl(":", coef_names)])
    ] %>% 
      pivot_longer(2:ncol(.), names_to = "term", values_to = "sd_x")
    # calculate linear combination for each climate
    effects <- bind_rows(
      biostat3::lincom(x, paste0(coef_names[1])),
      biostat3::lincom(x, paste0(c(coef_names[1], coef_names[2]), collapse = "+")),
      biostat3::lincom(x, paste0(c(coef_names[1], coef_names[3]), collapse = "+"))
    ) %>% 
      bind_rows() %>% rownames_to_column() %>% 
      separate(rowname, into = c("term", "climate"), sep = "\\+") %>% 
      mutate(
        climate = case_when(
          grepl("Hot", climate) ~ "Hot",
          grepl("Temperate", climate) ~ "Temperate",
          TRUE ~ "Cold"
        )
      )
    # standardise by climate
    effects_stan <- 
      effects %>% 
      merge(sd_y, by = "climate") %>% 
      merge(sd_x, by = c("climate", "term")) %>% 
      transmute(
        term, 
        Climate = factor(climate, levels = c("Hot", "Temperate", "Cold")),
        Estimate = Estimate * (sd_x / sd_y),
        conf.low = `2.5 %`  * (sd_x / sd_y),
        conf.high = `97.5 %`  * (sd_x / sd_y),
        statistic = Chisq, 
        p.value = `Pr(>Chisq)`
      )
    return(effects_stan)
  }) %>% 
  # combine different models into one data frame
  bind_rows() %>% 
  mutate(
    term = case_when(
      term %in% c("tmean_nf_12m") ~ "Temperature anomaly",
      term %in% c("tmean_nf_pos_12m") ~ "Temperature anomaly (+)",
      term %in% c("tmean_nf_neg_12m") ~ "Temperature anomaly (–)",
      term %in% c("tmean_xf_pos05_12m") ~ "Heat episode (temp.)",
      term %in% c("tmean_xf_neg05_12m") ~ "Cold episode (temp.)",
      term %in% c("utci_mean_xf_pos05_12m") ~ "Heat episode (UTCI)",
      term %in% c("utci_mean_xf_neg05_12m") ~ "Cold episode (UTCI)",
      term %in% c("spei3_neg_12m") ~ "Dry spell",
      term %in% c("spei3_pos_12m") ~ "Wet spell"
    ),
    dep_var = "Green vote share"
  )

# plots
figure3 <-
  lincom_concern %>% 
  bind_rows(lincom_vote) %>% 
  mutate(term = fct_rev(factor(term, levels = c(
    "Temperature anomaly", "Heat episode (temp.)", "Heat episode (UTCI)", "Dry spell"
  )))) %>% 
  ggplot(aes(x = Estimate, y = term)) +
  geom_pointrange(
    aes(xmin = conf.low, xmax = conf.high, colour = Climate, shape = Climate, fill = Climate), 
    position = position_dodge(0.5), size = 0.4
  ) +
  geom_vline(xintercept = 0, colour = "grey60") +
  scale_color_discrete(type = c(
    "Temperate" = "#d1bb8f",
    "Cold" = "#023a7a",
    "Hot" = "#7e3f00"
  )) +
  labs(x = "Standardized effect") +
  facet_wrap(vars(dep_var), ncol = 2) +
  theme(axis.title.y = element_blank())
ggsave("Figures/figure3.pdf", plot = figure3, width = 9, height = 5)

#### Figure S7: Dot-whisker plot of positive and negative anomalies by climate ####
xv <- c(
  "tmean_nf_pos_12m + tmean_nf_neg_12m +
  tmean_nf_pos_12m:climate + tmean_nf_neg_12m:climate",
  "tmean_xf_pos05_12m + tmean_xf_neg05_12m +
  tmean_xf_pos05_12m:climate + tmean_xf_neg05_12m:climate",
  "utci_mean_xf_pos05_12m + utci_mean_xf_neg05_12m +
  utci_mean_xf_pos05_12m:climate + utci_mean_xf_neg05_12m:climate",
  "spei3_pos_12m + spei3_neg_12m +
  spei3_pos_12m:climate + spei3_neg_12m:climate"
)
names_xv <- c(
  "Temperature anomaly (+)", "Temperature anomaly (–)",
  "Temperature anomaly (+) × Hot", "Temperature anomaly (+) × Temperate",
  "Temperature anomaly (–) × Hot", "Temperature anomaly (–) × Temperate",
  "Heat episode (temp.)", "Cold episode (temp.)",
  "Heat episode (temp.) × Hot", "Heat episode (temp.) × Temperate",
  "Cold episode (temp.) × Hot", "Cold episode (temp.) × Temperate",
  "Heat episode (UTCI)", "Cold episode (UTCI)",
  "Heat episode (UTCI) × Hot", "Heat episode (UTCI) × Temperate",
  "Cold episode (UTCI) × Hot", "Cold episode (UTCI) × Temperate",
  "Dry spell", "Wet spell",
  "Wet spell × Hot", "Wet spell × Temperate",
  "Dry spell × Hot", "Dry spell × Temperate"
)
tab_figure_s7 <- c(
  lapply(xv, function(i){ felm(
    f_felm(yvar = env_y, xvars = i, fevars = env_fv, clustervars = env_cv), 
    eb_climate, keepCX = TRUE, keepX = TRUE
  )}), 
  lapply(xv, function(i){ felm(
    f_felm(yvar = gv_y, xvars = i, fevars = gv_fv, clustervars = gv_cv),
    vote_climate, keepCX = TRUE, keepX = TRUE
  )})
)
tab_figure_s7_cs <- c(
  lapply(tab_figure_s7[1:length(xv)], conley_concern),
  lapply(tab_figure_s7[(1:length(xv))+length(xv)], conley_vote)
)

linc_concern <- 
  lapply(tab_figure_s7_cs[1:length(xv)], function(x){
    coef_names <- rownames(x$coefficients)
    # get SD of outcome by region 
    sd_y <- data.table(x$cY, climate = eb_climate$climate)[
      , .(sd_y = sd(env_concern)), by = .(climate)
    ]
    # get SD of covariates by region 
    sd_x <- data.table(x$cX, climate = eb_climate$climate)[
      , lapply(.SD, sd),
      by = .(climate),
      .SDcols = c(coef_names[!grepl(":", coef_names)])
    ] %>% 
      pivot_longer(2:ncol(.), names_to = "term", values_to = "sd_x")
    # calculate linear combination for each region
    effects <- bind_rows(
      biostat3::lincom(x, paste0(coef_names[1])),
      biostat3::lincom(x, paste0(c(coef_names[1], coef_names[3]), collapse = "+")),
      biostat3::lincom(x, paste0(c(coef_names[1], coef_names[4]), collapse = "+")),
      biostat3::lincom(x, paste0(coef_names[2])),
      biostat3::lincom(x, paste0(c(coef_names[2], coef_names[5]), collapse = "+")),
      biostat3::lincom(x, paste0(c(coef_names[2], coef_names[6]), collapse = "+"))
    ) %>% 
      bind_rows() %>% rownames_to_column() %>% 
      separate(rowname, into = c("term", "climate"), sep = "\\+") %>% 
      mutate(
        climate = case_when(
          grepl("Hot", climate) ~ "Hot",
          grepl("Temperate", climate) ~ "Temperate",
          TRUE ~ "Cold"
        )
      )
    # standardise by region
    effects_stan <- 
      effects %>% 
      merge(sd_y, by = "climate") %>% 
      merge(sd_x, by = c("climate", "term")) %>% 
      transmute(
        term, 
        Climate = factor(climate, levels = c("Hot", "Temperate", "Cold")),
        Estimate = Estimate * (sd_x / sd_y),
        conf.low = `2.5 %`  * (sd_x / sd_y),
        conf.high = `97.5 %`  * (sd_x / sd_y),
        statistic = Chisq, 
        p.value = `Pr(>Chisq)`
      )
    return(effects_stan)
  }) %>% 
  bind_rows() %>% 
  mutate(
    term = case_when(
      term %in% c("tmean_nf_pos_12m") ~ "Temperature anomaly (+)",
      term %in% c("tmean_nf_neg_12m") ~ "Temperature anomaly (–)",
      term %in% c("tmean_xf_pos05_12m") ~ "Heat episodes (temp.)",
      term %in% c("tmean_xf_neg05_12m") ~ "Cold episodes (temp.)",
      term %in% c("utci_mean_xf_pos05_12m") ~ "Heat episodes (UTCI)",
      term %in% c("utci_mean_xf_neg05_12m") ~ "Cold episodes (UTCI)",
      term %in% c("spei3_neg_12m") ~ "Dry spells",
      term %in% c("spei3_pos_12m") ~ "Wet spells"
    ),
    dep_var = "Environmental concern"
  )

linc_vote <- 
  lapply(tab_figure_s7_cs[(1:length(xv))+length(xv)], function(x){
    coef_names <- rownames(x$coefficients)
    # get SD of outcome by region 
    sd_y <- data.table(x$cY, climate = vote_climate$climate)[
      , .(sd_y = sd(green_vote)), by = .(climate)
    ]
    # get SD of covariates by region 
    sd_x <- data.table(x$cX, climate = vote_climate$climate)[
      , lapply(.SD, sd),
      by = .(climate),
      .SDcols = c(coef_names[!grepl(":", coef_names)])
    ] %>% 
      pivot_longer(2:ncol(.), names_to = "term", values_to = "sd_x")
    # calculate linear combination for each region
    effects <- bind_rows(
      biostat3::lincom(x, paste0(coef_names[1])),
      biostat3::lincom(x, paste0(c(coef_names[1], coef_names[3]), collapse = "+")),
      biostat3::lincom(x, paste0(c(coef_names[1], coef_names[4]), collapse = "+")),
      biostat3::lincom(x, paste0(coef_names[2])),
      biostat3::lincom(x, paste0(c(coef_names[2], coef_names[5]), collapse = "+")),
      biostat3::lincom(x, paste0(c(coef_names[2], coef_names[6]), collapse = "+"))
    ) %>% 
      bind_rows() %>% rownames_to_column() %>% 
      separate(rowname, into = c("term", "climate"), sep = "\\+") %>% 
      mutate(
        climate = case_when(
          grepl("Hot", climate) ~ "Hot",
          grepl("Temperate", climate) ~ "Temperate",
          TRUE ~ "Cold"
        )
      )
    # standardise by region
    effects_stan <- 
      effects %>% 
      merge(sd_y, by = "climate") %>% 
      merge(sd_x, by = c("climate", "term")) %>% 
      transmute(
        term, 
        Climate = factor(climate, levels = c("Hot", "Temperate", "Cold")),
        Estimate = Estimate * (sd_x / sd_y),
        conf.low = `2.5 %`  * (sd_x / sd_y),
        conf.high = `97.5 %`  * (sd_x / sd_y),
        statistic = Chisq, 
        p.value = `Pr(>Chisq)`
      )
    return(effects_stan)
  }) %>% 
  bind_rows() %>% 
  mutate(
    term = case_when(
      term %in% c("tmean_nf_pos_12m") ~ "Temperature anomaly (+)",
      term %in% c("tmean_nf_neg_12m") ~ "Temperature anomaly (–)",
      term %in% c("tmean_xf_pos05_12m") ~ "Heat episodes (temp.)",
      term %in% c("tmean_xf_neg05_12m") ~ "Cold episodes (temp.)",
      term %in% c("utci_mean_xf_pos05_12m") ~ "Heat episodes (UTCI)",
      term %in% c("utci_mean_xf_neg05_12m") ~ "Cold episodes (UTCI)",
      term %in% c("spei3_neg_12m") ~ "Dry spells",
      term %in% c("spei3_pos_12m") ~ "Wet spells"
    ),
    dep_var = "Green vote share"
  )

# plots
figure_s7 <-
  linc_concern %>% 
  bind_rows(linc_vote) %>% 
  mutate(
    term = fct_rev(factor(term, levels = c(
      "Temperature anomaly (+)", "Temperature anomaly (–)",
      "Heat episodes (temp.)", "Cold episodes (temp.)",
      "Heat episodes (UTCI)", "Cold episodes (UTCI)",
      "Dry spells", "Wet spells"
    )))
  ) %>% 
  ggplot(aes(x = Estimate, y = term)) +
  geom_pointrange(
    aes(xmin = conf.low, xmax = conf.high, colour = Climate, shape = Climate, fill = Climate), 
    position = position_dodge(0.5), size = 0.4
  ) +
  geom_vline(xintercept = 0, colour = "grey60") +
  scale_color_discrete(type = c(
    "Temperate" = "#d1bb8f",
    "Cold" = "#023a7a",
    "Hot" = "#7e3f00"
  )) +
  labs(x = "Standardized effect") +
  facet_wrap(vars(dep_var), ncol = 2) +
  theme(axis.title.y = element_blank())
# dw_region_concern
ggsave("Figures/figure_s7.pdf", plot = figure_s7, width = 9, height = 7)

#### Table S27: Climatic impacts on concerns and voting by income levels and changes ####
# standardize variables in concern dataset for interaction models
# see Table S5 for within_sd_concern
between_vars <- c(
  "gdp_ln_mean", "agri_frac", "irrigated_mean", "share_edu3_mean", "pop_share_u35", "LegislativeFractionalization_mean"
)
within_vars_concern <- c(
  "env_concern", "tmean_nf_12m", "gdp_ln_within", "CenterOfGravity_1"
)
eb_climate_stan <- copy(eb_climate)
eb_climate_stan[
    , c(between_vars, within_vars_concern) := lapply(.SD, function(x){(x - mean(x, na.rm = TRUE))}),
    .SDcols = c(between_vars, within_vars_concern)
  ][
    , c(between_vars) := lapply(.SD, function(x){(x / sd(x, na.rm = TRUE))}),
    .SDcols = c(between_vars)
  ][, `:=`(
    env_concern = env_concern / within_sd_concern[within_sd_concern$variable == "env_concern", "Within SD"],
    tmean_nf_12m = tmean_nf_12m / within_sd_concern[within_sd_concern$variable == "tmean_nf_12m", "Within SD"],
    gdp_ln_within = gdp_ln_within / within_sd_concern[within_sd_concern$variable == "gdp_ln_within", "Within SD"],
    CenterOfGravity_1 = CenterOfGravity_1 / within_sd_concern[within_sd_concern$variable == "CenterOfGravity_1", "Within SD"]
  )]

# standardize variables in voting dataset for interaction models
# see Table S6 for within_sd
within_vars_vote <- c(
  "green_vote", "tmean_nf_12m", "gdp_ln_within", "CenterOfGravity_1"
)
vote_climate_stan <- copy(vote_climate)
vote_climate_stan[
    , c(between_vars, within_vars_vote) := lapply(.SD, function(x){(x - mean(x, na.rm = TRUE))}),
    .SDcols = c(between_vars, within_vars_vote)
  ][
    , c(between_vars) := lapply(.SD, function(x){(x / sd(x, na.rm = TRUE))}),
    .SDcols = c(between_vars)
  ][, `:=`(
    green_vote = green_vote / within_sd_vote[within_sd_vote$variable == "green_vote", "Within SD"], 
    tmean_nf_12m = tmean_nf_12m / within_sd_vote[within_sd_vote$variable == "tmean_nf_12m", "Within SD"],
    gdp_ln_within = gdp_ln_within / within_sd_vote[within_sd_vote$variable == "gdp_ln_within", "Within SD"],
    CenterOfGravity_1 = CenterOfGravity_1 / within_sd_vote[within_sd_vote$variable == "CenterOfGravity_1", "Within SD"]
  )]

xv <- c(
  "tmean_nf_12m + gdp_ln_within + tmean_nf_12m:gdp_ln_within + tmean_nf_12m:gdp_ln_mean",
  "tmean_nf_12m + gdp_ln_within + tmean_nf_12m:gdp_ln_within + tmean_nf_12m:gdp_ln_mean + 
  tmean_nf_12m:agri_frac + tmean_nf_12m:irrigated_mean",#
  "tmean_nf_12m + gdp_ln_within + tmean_nf_12m:gdp_ln_within + tmean_nf_12m:gdp_ln_mean + 
  tmean_nf_12m:share_edu3_mean + tmean_nf_12m:pop_share_u35 + tmean_nf_12m:pop_urban",
  "tmean_nf_12m + gdp_ln_within + tmean_nf_12m:gdp_ln_within + tmean_nf_12m:gdp_ln_mean +
  tmean_nf_12m:LegislativeFractionalization_mean + CenterOfGravity_1 + tmean_nf_12m:CenterOfGravity_1",
  "tmean_nf_12m + gdp_ln_within + tmean_nf_12m:gdp_ln_within + tmean_nf_12m:gdp_ln_mean +
  tmean_nf_12m:climate"
)
names_xv <- c(
  "Temperature anomaly",
  "Log GDP per capita (within)",
  "Incumbent’s ideology (within)",
  "Temp. × GDP (within)",
  "Temp. × GDP (between)",
  "Temp. × agricultural fraction (between)",
  "Temp. × fraction irrigation (between)",
  "Temp. × share secondary education (between)",
  "Temp. × share under 35 (between)",
  "Temp. × urban region (between)",
  "Temp. × legislative fractionalization (between)",
  "Temp. × incumbent’s ideology",
  "Temp. × Hot climate (ref. cold)", 
  "Temp. × Temperate climate"
)
table_s27 <- c(
  lapply(xv, function(i){ felm(
    f_felm(yvar = env_y, xvars = i, fevars = env_fv, clustervars = env_cv), 
    eb_climate_stan, keepCX = TRUE, keepX = TRUE
  )}), 
  lapply(xv, function(i){ felm(
    f_felm(yvar = gv_y, xvars = i, fevars = gv_fv, clustervars = gv_cv),
    vote_climate_stan, keepCX = TRUE, keepX = TRUE
  )})
)
table_s27_cs <- c(
  lapply(table_s27[1:length(xv)], conley_concern),
  lapply(table_s27[(1:length(xv))+length(xv)], conley_vote)
)
export_html(table_s27_cs, out = "Tables/table_s27.html", cov_labels = names_xv)#

#### Figure 3: Effects of temperature anomaly on by income level and climate ####
# available on request

#### Table S28: Effects of environmental concern on Green voting ####
xv <- c(
  "Env. concern (12 months lagging mean)" = "env_concern_12m",
  "Env. concern (24 months lagging mean)" = "env_concern_24m",
  "Env. concern (12 months leading mean)" = "env_concern_12m_lead",
  "Env. concern (24 months leading mean)" = "env_concern_24m_lead"
)
fit <- c("Env. concern (predicted values, 12 month lagging mean)" = "env_concern_12m")
iv <- c(
  "Temperature anomaly" = "tmean_nf_12m",
  "Heat episode (temp.)" = "tmean_xf_pos05_12m",
  "Heat episode (UTCI)" = "utci_mean_xf_pos05_12m",
  "Dry spell" = "spei3_neg_12m"
)
table_s28 <- c(
  lapply(xv, function(i){ felm(
    f_felm(yvar = gv_y, xvars = i, fevars = gv_fv, clustervars = gv_cv),
    vote_climate, keepCX = TRUE, keepX = TRUE,  exactDOF = TRUE
  )}),
  lapply(iv, function(i){ felm(
    f_felm(yvar = gv_y, fit = fit, ivars = i, fevars = gv_fv, clustervars = gv_cv),
    vote_climate, keepCX = TRUE, keepX = TRUE,  exactDOF = TRUE
  )})
)
table_s28_cs <- c(
  lapply(table_s28, conley_vote)
)
export_html(
  lapply(table_s28_cs, standardise), 
  out = "Tables/table_s28_stage2.html",
  cov_labels = c(names(xv), names(fit)),
  dep_labels = "Green vote share",
  add_lines = list(
    c("F statistic", rep("", length(xv)), sapply(table_s28[5:(4+length(iv))], get_f_iv)),
    c("Unit fixed effects", rep("×", length(xv) + length(iv))),
    c("Period fixed effects", rep("×", length(xv) + length(iv))),
    c("Spatial cutoff (km)", rep(gv_spc, length(xv) + length(iv))),
    c("Temporal cutoff (years)", rep(gv_tc, length(xv) + length(iv))),
    c("Observations", sapply(table_s28_cs, get_n)),
    c("Overall R²", sapply(table_s28_cs[1:4], get_r2_overall)),
    c("Within R²", sapply(table_s28_cs[1:4], get_r2_within))
  )
)
quietly(stargazer(
  lapply(table_s28_cs[5:(4+length(iv))], standardise, stage1 = TRUE),
  out = "Tables/table_s28_stage1.html", type = "html", covariate.labels = names(iv),
  dep.var.labels = "Environmental concern",
  add.lines = list(
    c("Within R²", sapply(table_s28_cs[5:(4+length(iv))], get_r2_within, stage1 = TRUE))
  )
))

#### Table S29: Models testing for the impacts of changes in the concern measurement on response behavior ####
# available on request

#### Table S30: Matrix of bivariate correlations between climate variables ####
cor_vars <- c(
  "tmean_nf_12m" = "Temperature anomaly",
  "tmean_nf_pos_12m" = "Temperature anomaly (positive)",
  "tmean_nf_neg_12m" = "Temperature anomaly (negative)",
  "tmean_xf_pos05_12m" = "Heat episode (temp.)",
  "tmean_xf_neg05_12m" = "Cold episode (temp.)",
  "utci_mean_xf_pos05_12m" = "Heat episode (UTCI)",
  "utci_mean_xf_neg05_12m" = "Cold episode (UTCI)",
  "spei3_neg_12m" = "Dry spell (SPEI3)",
  "spei3_pos_12m" = "Wet spell (SPEI3)",
  "flood_frac_12m" = "Flooding",
  "burn_frac_12m" = "Wildfire"
)
cor_vars_names <- names(cor_vars)
table_s30 <- cor(eb_climate[, ..cor_vars_names], use = "pair")
rownames(table_s30) <- cor_vars
colnames(table_s30) <- cor_vars
print(
  xtable(table_s30),
  file = "Tables/table_s30.html",
  type = "html"
)

